function [R2,alphas,betas] = solveR2(alphas0,betas0,w,tspan,dt)
% 
tab_w = [0.10,0.15,0.20,0.50,1.00];
tab_nu = [0.0451e-4,1.3915e-005,0.3583e-4,0.1232e-3,0.3939e-3];
tab_sigma2 = [1.0810e-3,6.2341e-3,2.1983e-2,0.0990,0.3007];

% 
% assert((w>=min(tab_w)&& w<=max(tab_w)),'illegal w = %f\n',w);
assert((w>=0 && w<=max(tab_w)),'illegal w = %f\n',w);
    
nu = interp1(tab_w,tab_nu,w,'linear','extrap');
sigma2 = interp1(tab_w,tab_sigma2,w,'linear','extrap');
sigma = sqrt(sigma2);

alphas = alphas0;
betas = betas0;

for n = 1:length(alphas)
    lambda  = -nu*n^2;
    alphas(n) = solveAlpha(alphas(n),lambda,sigma,tspan,dt);
    betas(n) = solveAlpha(betas(n),lambda,sigma,tspan,dt);
end
R2 = sum(alphas.^2+betas.^2);
